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ABSTRACT 

The zonal flow in Jupiter's upper troposphere is organized into alternating retrograde and prograde jets, with a prograde (super- 
rotating) jet at the equator. Existing models posit as the driver of the flow either differential radiative heating of the atmosphere or 
intrinsic heat fluxes emanating from the deep interior; however, they do not reproduce all large-scale features of Jupiter's jets and 
thermal structure. Here it is shown that the difficulties in accounting for Jupiter's jets and thermal structure resolve if the effects 
of differential radiative heating and intrinsic heat fluxes are considered together, and if upper-tropospheric dynamics are linked to a 
magnetohydrodynamic (MHD) drag that acts deep in the atmosphere and affects the zonal flow away from but not near the equator. 
Baroclinic eddies generated by differential radiative heating can account for the off-equatorial jets; meridionally propagating equatorial 
Rossby waves generated by intrinsic convective heat fluxes can account for the equatorial superrotation. The zonal flow extends deeply 
into the atmosphere, with its speed changing with depth, away from the equator up to depths at which the MHD drag acts. The theory 
is supported by simulations with an energetically consistent general circulation model of Jupiter's outer atmosphere. A simulation that 
incorporates differential radiative heating and intrinsic heat fluxes reproduces Jupiter's observed jets and thermal structure and makes 
testable predictions about as-yet unobserved aspects thereof. A control simulation that incorporates only differential radiative heating 
but not intrinsic heat fluxes produces off-equatorial jets but no equatorial superrotation; another control simulation that incorporates 
only intrinsic heat fluxes but not differential radiative heating produces equatorial superrotation but no off-equatorial jets. The proposed 
mechanisms for the formation of jets and equatorial superrotation likely act in the atmospheres of all giant planets. 



1. Introduction 

The zonal flow in Jupiter's upper troposphere has been 
inferred by tracking cloud features, which move with 
the horizontal flow in the layer between about 0.5 and 
1 bar atmospheric pressure (Ingersoll et al. 2004; West 
et al. 2004; Vasavada and Showman 2005). In this layer, 
the zonal flow is organized into a strong prograde (su- 
perrotating) equatorial jet and an alternating sequence 
of retrograde and prograde off-equatorial jets (Fig. la). 
This flow pattern has been stable at least between the 
observations by the Voyager and Cassini spacecrafts in 
1979 and 2000, with some variations in jet speeds, for 
example, a slowing of the prograde jet at 21°N plane- 
tocentric latitude by ~40ms _1 (Porco et al. 2003; In- 
gersoll et al. 2004). The zonal flow in layers above the 
clouds has been inferred from the thermal structure of 
the atmosphere, using the thermal wind relation between 
meridional temperature gradients and vertical shears of 
the zonal flow. Meridional temperature gradients and 
thus vertical shears between 0.1 and 0.5 bar are gener- 
ally small: meridional temperature contrasts along iso- 
bars do not exceed ^10 K (Conrath et al. 1998; Simon- 
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Miller et al. 2006). The thermal stratification in the same 
layer is statically stable (Simon-Miller et al. 2006; Read 
et al. 2006). About the zonal flow in lower layers, it is 
only known that at one site at 6.4°N planetocentric lat- 
itude, where the Galileo probe descended into Jupiter's 
atmosphere, it is prograde and increases with depth from 
^90ms _1 at 0.7 bar to ^170ms _1 at 4 bar; beneath, 
it is relatively constant up to at least ~20 bar (Atkin- 
son et al. 1998). At the same site, the thermal stratifi- 
cation is statically stable but approaches neutrality with 
increasing depth between 0.5 and 1.7 bar; beneath, it is 
statically nearly neutral or neutral up to at least ^20 bar 
(Magalhaes et al. 2002). These are the large-scale fea- 
tures (if the Galileo probe data are representative of large 
scales) that a minimal model of Jupiter's general circu- 
lation should be able to reproduce. 

There are two plausible energy sources for Jupiter's 
general circulation. First, ~8Wm -2 of solar radia- 
tion are absorbed in Jupiter's atmosphere (Hanel et al. 
1981), where the time-mean insolation at the top of 
the atmosphere, given Jupiter's small obliquity of 3°, 
varies approximately with the cosine of latitude. Sec- 
ond, ~6Wm -2 of intrinsic heat fluxes emanate from 
Jupiter's deep interior (Ingersoll et al. 2004; Guillot et al. 
2004; Guillot 2005); observations of convective storms 
(Gierasch et al. 2000; Porco et al. 2003; Sanchez -Lavega 
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Zonal velocity (m s~ 1 ) 

FIGURE 1: Zonal flow in Jupiter's upper troposphere and in 
simulations, (a) Zonal velocity on Jupiter, inferred by tracking 
cloud features from the Cassini spacecraft (Porco et al. 2003) 
(orange) and in Jupiter simulation at 0.65 bar (blue), (b) Zonal 
velocity at 0.65 bar in control simulations: with intrinsic heat 
fluxes but with uniform insolation at the top of the atmosphere 
(magenta), and with differential insolation but without intrin- 
sic heat fluxes (light blue). Zonal velocities in simulations are 
zonal and temporal means in statistically steady states (over 
1500 days for Jupiter simulation and over 900 days for control 
simulations); differences between the (statistically identical) 
hemispheres here and in subsequent figures are indicative of 
sampling variability. For Jupiter, latitude here and throughout 
this paper is plane tocentric; the simulated planets are spherical, 
so planetocentric and planetographic latitudes are identical. 

et al. 2008) and the neutral or nearly neutral thermal 
stratification along much of the Galileo probe descent 
path show that the intrinsic heat fluxes are at least par- 
tially convective. Existing models of Jupiter's zonal flow 
posit as the driver of the general circulation either dif- 
ferential radiative heating of the atmosphere or intrinsic 
convective heat fluxes (e.g., Busse 1976, 1994; Williams 
1979; Ingersoll et al. 2004; Vasavada and Showman 
2005); however, they do not reproduce all large-scale 
features of the jets and thermal structure. For exam- 
ple, in parameter regimes relevant for Jupiter, they gen- 
erally do not produce equatorial superrotation, unless 
artifices are employed such as imposing an additional 
heat source near the equator (Williams 2003) or assum- 
ing excessively viscous flow and permitting intrinsic heat 
fluxes several orders of magnitude stronger than Jupiter's 
(Heimpel et al. 2005; Heimpel and Aurnou 2007). 



Here we show that the difficulties in accounting for 
Jupiter's jets and thermal structure resolve if the effects 
of differential radiative heating and intrinsic heat fluxes 
are considered together, and if upper-tropospheric dy- 
namics are linked to drag that acts deep in the atmo- 
sphere and affects the zonal flow away from but not near 
the equator. The key is to distinguish the different ways 
in which eddies can be generated near the equator and 
away from it, and to consider their role and the role of 
drag at depth in the balance of angular momentum (the 
angular momentum component about the planet's spin 
axis). 

First we describe how eddies can be generated near 
the equator and away from it, how convectively gener- 
ated equatorial waves can lead to equatorial superrota- 
tion, and how the angular momentum balance of the up- 
per troposphere is linked to drag and constrains the flow 
at depth (sections 2-A). Then we use simulations with a 
three-dimensional general circulation model (GCM) of 
a thin shell in Jupiter's outer atmosphere, with an ide- 
alized representation of effects of drag deep in the at- 
mosphere, to demonstrate that the mechanisms proposed 
can account for Jupiter's observed jets and thermal struc- 
ture (sections 5-7). 

2. Eddy generation and angular momentum fluxes 

Tracking of cloud features shows that eddies in Jupiter's 
upper troposphere transport angular momentum merid- 
ionally from retrograde into prograde jets; the eddy an- 
gular momentum fluxes extend at least over the layer be- 
tween about 0.5 and 1 bar, and imply there a mean con- 
version rate from eddy to mean flow kinetic energy per 
unit mass of order 10~ 5 to 10~ 4 W kg -1 (Ingersoll et al. 
1981; Salyk et al. 2006). If the eddy angular momen- 
tum fluxes extended unabatedly over a layer more than 
tens of bars thick, the total conversion rate from eddy 
to mean flow kinetic energy would exceed the rate at 
which the atmosphere takes up energy from solar radi- 
ation and intrinsic heat fluxes combined. For example, 
if the eddy angular momentum fluxes extended over a 
layer of 2.5 bar thickness, and if vertical variations of 
the zonal flow over this layer are negligible, the total 
conversion rate would already be of order 1 Wm~ 2 — 
that is, it would amount to ^5% of the energy uptake by 
the atmosphere (Salyk et al. 2006). But probably only 
a small fraction of the atmosphere's energy uptake is 
available to generate eddy kinetic energy. This means 
that columnar-flow models of Jupiter's general circula- 
tion, in which the momentum exchange between eddies 
and the mean flow extends over layers with thickness of 
order 10 6 bar and greater, are not viable on energetic 
grounds: they require energy sources several orders of 
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magnitude larger than Jupiter's (e.g., Sun et al. 1993; 
Heimpel et al. 2005; Heimpel and Aurnou 2007). It also 
means that the thin-shell approximation — in which the 
distance from any point in the atmosphere to the planet's 
center is taken to be constant and equal to the planetary 
radius — is adequate for and will be made in the follow- 
ing considerations of the tropospheric eddy transport of 
angular momentum (although the zonal flow can extend 
deeply, see section 4). 

Meridional eddy transport of angular momentum 
is evidence of meridional eddy propagation and irre- 
versible dynamics (e.g., Edmon et al. 1980). In a thin 
shell, eddies that propagate meridionally transport angu- 
lar momentum in the opposite direction of their prop- 
agation if the meridional gradient of absolute vorticity, 
or depth-averaged potential vorticity, is positive (north- 
ward) between their generation and dissipation (break- 
ing) regions. In the dissipation regions, irreversible 
meridional mixing of absolute vorticity then leads to 
southward eddy vorticity fluxes and divergence of merid- 
ional eddy angular momentum fluxes; the compensat- 
ing convergence of meridional eddy angular momentum 
fluxes occurs in the generation regions, implying angu- 
lar momentum transport from the dissipation into the 
generation regions (Kuo 1951; Held 1975; Andrews and 
Mclntyre 1976, 1978; Rhines 1994; Vallis 2006, chap- 
ter 12). In Jupiter's upper troposphere, the meridional 
gradient of absolute vorticity (and potential vorticity) is 
generally positive, except in narrow latitude bands at the 
centers of retrograde jets (Ingersoll et al. 1981; Read 
et al. 2006). The direction of the observed eddy angu- 
lar momentum fluxes hence indicates generation of ed- 
dies in prograde jets, dissipation in retrograde jets, and 
propagation in between. 

In prograde off-equatorial jets, eddies can be gener- 
ated by baroclinic instability. Prograde jets are baroclin- 
ically more unstable than retrograde jets if the speed (ab- 
solute velocity value) of the zonal flow in prograde jets 
decreases with depth and in retrograde jets either also 
decreases with depth or, with weaker vertical shear, in- 
creases with depth. By thermal wind balance, merid- 
ional temperature gradients along isobars then are equa- 
torward in prograde jets and either are also equatorward 
but weaker or are poleward in retrograde jets. Indeed, in 
Jupiter's upper troposphere (at ^0.25 bar), meridional 
temperature gradients have been observed to be equator- 
ward in prograde off-equatorial jets and either equator- 
ward but weaker or poleward in retrograde off-equatorial 
jets [compare Fig. 6 in Simon-Miller et al. (2006) with 
Fig. 3 in Vasavada and Showman (2005)]. Thus, merid- 
ional propagation of eddies generated by baroclinic in- 
stability preferentially in prograde jets can account for 
the angular momentum transport from retrograde into 



prograde off-equatorial jets, as demonstrated in baro- 
clinic GCMs that produce multiple jets (e.g., Williams 
1979; O'Gorman and Schneider 2008). 

In the prograde equatorial jet, eddies can be gener- 
ated by intrinsic convective heat fluxes. Unlike in higher 
latitudes, Coriolis forces are small in the equatorial re- 
gion, where the Rossby number is order one or greater. 
Horizontal pressure gradients there are limited by the 
necessity to be balanced primarily by inertial accelera- 
tions, rather than or in addition to Coriolis accelerations. 
Hence, they are small (of order Froude number), and 
so are horizontal temperature gradients on scales large 
enough that hydrostatic balance holds (Charney 1963). 1 
On such large scales, therefore, diabatic heating Q can- 
not be balanced by temperature fluctuations or horizontal 
temperature advection, as in higher latitudes. Instead, as 
in the tropics of Earth's atmosphere, it is primarily bal- 
anced by the adiabatic cooling associated with vertical 
motion u acting on a (possibly small) static stability S, 
giving the weak temperature gradient approximation of 
the thermodynamic equation 2 (Held and Hoskins 1985; 
Sobel et al. 2001) 

uS^-Q. (1) 

By mass conservation, • v x = —d p uj, vertical gradi- 
ents of diabatic heating and/or static stability thus induce 
large-scale horizontal divergence 

V h -v x *d p (Q/S), (2) 

where we have decomposed the horizontal velocity v = 
v x + v$ into divergent (v x ) and rotational (v$) com- 

1 The angular momentum and hydrostatic equations imply that, near 
the equator, horizontal variations in pressure p and potential tempera- 
ture 9 scale as 5p/p ~ 59/9 ~ Fr, where Fr = UV/{gH) is a 
Froude number, H is the scale height, U is a mean zonal velocity 
scale, and V is the greater of an eddy velocity scale or mean merid- 
ional velocity scale (cf. Charney 1963). For Jupiter parameters and 
with H ~ 20km, U ~ 100ms- 1 , and V ~ 10ms- 1 (Porco et al. 
2003; Salyk et al. 2006), one obtains Fr ~ 10~ 3 . This scaling holds 
where the Rossby number Ro = U/\fL\, with length scale of flow 
variations L, is order one or greater. With f = /3y and L < \y\, 
it follows that the scaling holds at least within meridional distances 
\y\ ~ (C///3) 1 / 2 of the equator (Sobel et al. 2001), which, for Jupiter, 
is within ~5000 km or ~4° of the equator. 

2 We use pressure coordinates, with S = —d p 9 and ui = Dp/ Dt. 
To obtain the relative magnitudes of the terms on the left-hand side of 
the thermodynamic equation (dt+v ■ Vj,)8-tjS = Q, withisobaric 
horizontal derivative operator V/j, we assume the explicit time deriva- 
tive and horizontal advection terms scale as (V/L)S9, and the vertical 
advection term scales as (V/L)A9. We have used mass conservation 
to relate vertical to horizontal velocity scales, and A9 is a vertical po- 
tential temperature change over a scale height. The vertical advection 
term then dominates, and the weak temperature gradient approxima- 
tion (1) is adequate, if vertical potential temperature changes satisfy 
A9/9 3> 59/9 ~ Fr ~ 10~ 3 . In Jupiter's equatorial troposphere, 
if one takes the thermal stratification along the Galileo probe descent 
path as representative, this is assured at least above 1.7 bar (Magalhaes 
et al. 2002). 
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ponents. As discussed by Sardeshmukh and Hoskins 
(1988), this horizontal divergence drives a rotational 
flow: Vortex stretching and vorticity advection by the 
divergent flow are a source 

R = -Ca(V h ■ V X ) - (V X • V h )Ca 0) 

of rotational flow, as can be seen from the equation for 
the absolute vorticity £ a = / + k • (V^ x v$) = / + C 
in the equatorial troposphere in the form 

+ v* ■ V h )C. « (4) 

Neglected here are the baroclinic term, consistent with 
the weak temperature gradient approximation (Charney 
1963; Sobel et al. 2001), as well as friction and the verti- 
cal advection and tilting terms. It follows that convective 
heating fluctuations that cannot be balanced by slow ra- 
diative processes induce fluctuations in the large-scale 
horizontal divergence, and these represent a source 

Rl = R-R= -V fc • (Cav x - (a^x) (5) 

of vorticity fluctuations and thus a source of Rossby 
waves (Sardeshmukh and Hoskins 1988). (Overbars de- 
note isobaric zonal and temporal means and primes de- 
viations therefrom.) 

Horizontally, the Rossby wave source (5) can be ex- 
pected to be large in the equatorial region because the 
Rossby number there is order one or greater, and large- 
scale horizontal flow fluctuations induced by convective 
heating fluctuations are divergent at leading order. In 
contrast, the Rossby number in higher latitudes is small, 
large-scale horizontal flow fluctuations are nondivergent 
at leading order, and baroclinic generation of vortic- 
ity fluctuations, neglected in the Rossby wave source 
(5), can be important. Scale analysis 3 suggests that the 
Rossby wave source (5) is largest in the prograde equa- 
torial jet and in the strong prograde jet at 21°N — that 

3 For the scale analysis, we use the eddy velocity scale V ~ 
10 ms and the length scale of flow variations L ~ 5000 km both 
for eddies and for mean fields (cf. footnote 1); we take the scales to 
be invariant with latitude, as indicated, for the velocity scale, by the 
Cassini image analysis of Salyk et al. (2006). Then, where Ro > 1 
and if horizontal velocity fluctuations are divergent at leading order, 
the Rossby wave source owing to advection of planetary vorticity by 
the divergent flow is of order R' ~ /3V ~ 5 X 10~ n s~ 2 ; the Rossby 
wave source owing to stretching of planetary vorticity, evaluated at 4° 
latitude, is of the same order R' ~ \ fV\/L ~ 5x lO" 11 s~ 2 . Where 
Ro < 1, the divergent velocity is of order Ro V, and the Rossby wave 
source is of order R' ~ Ro \ fV\/L = UV/L 2 . This is of the same 
order (R' ~ 5 X 10~ n s~ 2 ) as the Rossby wave source near the equa- 
tor where the mean zonal velocity scale is of order U ~ 100 ms -1 
(throughout the prograde equatorial jet and in the strong prograde jet 
at 21°N, see Fig. la). At other latitudes, the mean zonal velocity scale 
(U < 30 ms -1 ) and the Rossby wave source (R' < 10~ n s~ 2 ) are 
smaller, albeit only by O(l) factors. 



is, not only in the latitude band in which the Rossby 
number is order one or greater (within ^4° of the equa- 
tor, see footnote 1), but in a slightly wider latitude band 
around the equator and in strong jets generally. Verti- 
cally, the Rossby wave source (5) should be largest in 
the upper troposphere, below the top of the convective 
outflows, where fluctuations in the large-scale horizon- 
tal divergence (2) can be expected to be largest: Con- 
vective heating fluctuations and their vertical gradients 
can be expected to be largest there, and the static sta- 
bility has substantial vertical gradients (Magalhaes et al. 
2002), marking the transition from lower layers that are 
neutrally stratified by convection to upper layers that are 
more stably stratified because the stabilizing radiative 
heating from above begins to have an effect. The Rossby 
wave source (5), then, should be largest in Jupiter's equa- 
torial upper troposphere if convective heating fluctua- 
tions there are sufficiently strong. 

3. Generation of equatorial superrotation 

By Hide's theorem, the eddy transport of angular mo- 
mentum into the equatorial region observed in Jupiter's 
upper troposphere is necessary for the existence of a pro- 
grade equatorial jet (Hide 1969; Schneider 1977; Sara- 
vanan 1993; Held 1999; Schneider 2006). To see how 
such angular momentum transport can come about, con- 
sider the eddy enstrophy equation implied by the vortic- 
ity equation (4), 

d t (^/2+^cd v ( a ^we, (6) 

where y — a<p is the meridional coordinate with plan- 
etary radius a and latitude <f>. Neglected here is the ad- 
vection of eddy enstrophy by the eddies themselves, a 
triple correlation term that seems to be about an order 
of magnitude smaller in Jupiter's equatorial upper tro- 
posphere than the retained term v'^C d y ( a . 4 In a sta- 
tistically steady state, upon division by d y Q a ^ and 
multiplication by the thin-shell approximation of the dis- 
tance r± — acoscj) to the planet's spin axis, the eddy 
enstrophy equation becomes the wave activity balance 
(Andrews and Mclntyre 1976, 1978; Edmon et al. 1980) 



Gwt4C / r± = -div(uV^r ± ). (7) 

4 To estimate the relative magnitude of the terms, we roughly ap- 
proximate the absolute vorticity gradient d y £ a by (i (Read et al. 2006) 
and use the horizontal length scale of flow variations L ~ 5000 km 
(see footnote 3) and the rotational meridional eddy velocity scale V ~ 
10 ms" 1 (Salyk et al. 2006); that is, we assume the rotational and 
divergent velocities are of the same order (cf. footnote 3). Then, the 
magnitude of the neglected eddy enstrophy advection div(n^,f' 2 )/2 
relative to the retained term v' 9 C d y Ca is of order V/(f3L 2 ) ~ 0.1. 
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Here, div(-) is the isobaric meridional divergence oper- 
ator, and 



R'C 

OyCa 

represents the generation of wave activity A, with 



A = 



l C' 2 



2 dyCa 



rj_. 



(8) 



(9) 



The wave activity balance (7) states that at latitudes at 
which wave activity is generated (G > 0), the eddy 
vorticity flux v'^C' is directed northward, implying con- 
vergence of (rotational) eddy angular momentum fluxes 
u'yv'q, r± (per unit mass); conversely, at latitudes at 
which wave activity is dissipated (G < 0), there is di- 
vergence of eddy angular momentum fluxes. The merid- 
ional eddy angular momentum flux is equal to minus 
the meridional wave activity flux. Thus, generation of 
wave activity in the equatorial region and radiation to 
and dissipation in higher latitudes entail angular momen- 
tum transport from higher latitudes into the equatorial 
region. See Andrews and Mclntyre (1976, 1978), Plumb 
(1979), Mclntyre (1980), Edmon et al. (1980), and Val- 
lis (2006, chapter 7) for further discussion and general- 
izations of these results from wave-mean flow interac- 
tion theory; and see Suarez and Duffy (1992) and Sara- 
vanan (1993) for numerical demonstrations that (station- 
ary) Rossby wave sources in the equatorial region entail 
convergence of eddy angular momentum fluxes and can 
generate superrotation. 

The wave activity balance (7) implies that in Jupiter's 
upper troposphere, convergence of eddy angular mo- 
mentum fluxes in the equatorial region is to be expected 
if Rossby wave generation by convective heating fluctu- 
ations is sufficiently strong. Contributions to the Rossby 
wave source R' that do not depend on vorticity fluctua- 
tions to the extent that d t (' ~ R' , can be expected 
to contribute positively to the eddy enstrophy generation 
R'C,' and thus, because the absolute vorticity gradient is 
positive (Read et al. 2006), to the wave activity genera- 
tion G. Some of the wave activity so generated may dis- 
sipate near or in its generation region, as contributions 
to the Rossby wave source R' that depend on vorticity 
fluctuations £' may damp them. The principal damping 
term is — (div v x )C' 2 ', as pointed out Sardeshmukh and 
Hoskins (1988), this term contributes negatively to the 
eddy enstrophy and wave activity generation where the 
mean meridional flow is divergent, for example, in the 
horizontal outflows of ascending branches of any (Eule- 
rian) mean meridional circulation cells. Near the equa- 
tor, however, the damping limits itself as the strength of 
convective heating fluctuations increases: Wave activ- 
ity generation near the equator entails convergence of 



eddy angular momentum fluxes, which implies weak- 
ened poleward or even equatorward mean meridional 
flow, as can be seen from the zonally and temporally 
averaged angular momentum balance in a statistically 
steady state, 



-u • VMn = fv x r± = S e + S n 



(10) 



Here, u is the three-dimensional velocity vector, Mq = 
Qr± is the angular momentum per unit mass owing to 
the planetary rotation (with constant angular velocity ft), 
and the first equality holds in the thin-shell approxima- 
tion. We have again neglected friction, we have used 
the fact that the mean meridional flow is irrotational 
(v = v x ), and 



S e = div(uV rj_) + d p (u'u)' r±), 



(11a) 



and 



div(uvr±) + d p (ucor±), (lib) 



are the divergences of fluxes of relative angular momen- 
tum per unit mass M u = ur± owing to eddies and mean 
meridional circulations. The angular momentum balance 
(10) shows that as the strength of convective heating fluc- 
tuations increases, increasing wave activity generation 
and the increasing eddy angular momentum flux conver- 
gence (decreasing S e ) it entails either weaken any di- 
vergence of the mean meridional flow at the equator or 
even lead to convergence and thus to eddy enstrophy and 
wave activity amplification. 5 Sufficiently strong convec- 
tive heating fluctuations, then, can be expected to lead to 
net wave activity generation (G > 0) and convergence of 
eddy angular momentum fluxes in the equatorial region. 

Convergence of eddy angular momentum fluxes in the 
equatorial region accelerates the zonal flow. As a con- 
sequence, a prograde equatorial jet forms if other pro- 
cesses that may decelerate it — e.g., drag at depth linked 
to upper-tropospheric dynamics through mean merid- 
ional circulations (see section 4) — are sufficiently weak. 
Such an equatorial jet can be expected to occupy a lati- 
tude band at least as wide as the meridional decay scale 
of long equatorial Rossby waves: the equatorial Rossby 
radius L e = (c/f3) 1 ^ 2 with tropospheric gravity wave 
speed c (Matsuno 1966; Gill 1982, chapter 11). For 
Jupiter, if one takes the speed c w 450 ms -1 of the 
waves observed after the impact of comet Shoemaker- 
Levy 9 as the relevant gravity wave speed (Ingersoll and 
Kanamori 1995) — a wave speed roughly consistent with 
the tropospheric thermal stratification along the Galileo 



5 This is one of several eddy-mean flow feedbacks that can lead to 
rapid transitions to superrotation as parameters such as the strength of 
the convective heating are varied. Other feedbacks are discussed by 
Saravanan (1993) and Held (1999). 
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probe descent path (cf. Magalhaes et al. 2002) — one ob- 
tains L e 9500 km, or a lower bound on the half-width 
of the jet of ^8°. This is similar to, albeit smaller than, 
the half-width of Jupiter's equatorial jet in the upper tro- 
posphere (Fig. la). 6 

An upper bound on the speed of the equatorial jet can 
be obtained by vorticity homogenization arguments. The 
jet speed is determined by the change in vorticity within 
the equatorial jet owing to the meridional redistribution 
of absolute vorticity relative to a state of solid body ro- 
tation. If the zonal flow is zero at a distance L s from 
the equator and if this distance is small enough that the 
small-angle approximation for latitudes is adequate, the 
jet speed at the equator is U w A(L S , where A£ is the 
absolute value of the relative vorticity averaged between 
the equator and a distance L s away from it. The absolute 
vorticity redistribution reaches an end state, and merid- 
ional Rossby wave propagation ceases, when the abso- 
lute vorticity within the jet is homogenized in each hemi- 
sphere. If the absolute vorticity at the zeros of the jet is 
approximated by ±/3L s , the state of homogenized ab- 
solute vorticity corresponding to maximum pro grade jet 
speed for a given jet width has absolute vorticity +(3L S 
in the northern and —(3L S in the southern hemisphere, 
with a barotropically stable jump at the equator: an abso- 
lute vorticity "staircase" (Mclntyre 1982; Dritschel and 
Mclntyre 2008). The absolute value of the relative vor- 
ticity averaged between the equator and a distance L s 
away from it then is A( w 0L s /2, and one obtains the 
upper bound on the jet speed 



If one further substitutes L s ~ L e , one obtains U < c/2, 
that is, a bound on the jet speed that only depends on the 
gravity wave speed. For Jupiter, with c w 450 ms -1 , 
this gives U < 225 ms -1 , which is of the same order 
of, albeit larger than, the observed equatorial jet speed in 
the upper troposphere (Fig. la); it is closer to the equa- 
torial jet speed at greater depth along the Galileo probe 
descent path. If the actual half-width of Jupiter's equa- 
torial jet is used (L s > L e ), the overestimation of the 
jet speed in the upper troposphere is more pronounced; 
a state of homogenized absolute vorticity in each hemi- 
sphere is not attained in the upper troposphere (Read 
et al. 2006). Nonetheless, these arguments show that 

6 Other lower bounds on the half-width of the equatorial jet are the 
half-widths of the latitude bands with (i) substantial fluctuations in the 
horizontal divergence or (ii) relative angular momentum transport by 
mean meridional circulations. Both latitude bands are characterized 
by Ro > 1 (e.g., Sobel et al. 2001; Schneider 2006). With the upper 
bound U < c/2 on the equatorial jet speed given below, one finds that 
the half- width of the latitude band with Ro > 1 is constrained to be of 
the same order of but smaller than L e (cf. footnote 1). 



the speed of prograde equatorial jets, to the extent that 
a state of homogenized absolute vorticity in each hemi- 
sphere is being approached, can be expected to increase 
roughly linearly with (3 and quadratically with their half- 
width. The speed of such jets does not depend directly, 
as is sometimes surmised, on the total energy uptake by 
the atmosphere or the total kinetic energy dissipation, but 
on emergent properties such as the thermal stratification, 
which determines the gravity wave speed. Consistent 
with these arguments and the similar radii and rotation 
rates of Jupiter and Saturn, Saturn's prograde equatorial 
jet is about 2 times wider and 3—4 times stronger than 
Jupiter's. This may be a consequence of a greater grav- 
ity wave speed on Saturn. 

4. Drag at depth and mean meridional circulations 

In a statistically steady state, eddy angular momentum 
transport from retrograde into prograde jets in the up- 
per troposphere must be balanced by angular momen- 
tum transport by mean meridional circulations or eddies 
in other layers or by drag on the zonal flow. Drag acts on 
the zonal flow deep in Jupiter's atmosphere, where hy- 
drogen, its primary constituent, undergoes a continuous 
transition from an electrically semi-conducting molec- 
ular phase in the outer layers to a conducting metallic 
phase in the interior (Guillot et al. 2004; Guillot 2005). 
Where the atmosphere is electrically conducting, any 
flow advects the magnetic field and induces an electric 
current. The Ohmic dissipation of the induced current 
implies a magnetohydrodynamic (MHD) drag on the 
flow (Grote and Busse 2001; Liu et al. 2008). This MHD 
drag we take to provide the principal angular momentum 
dissipation mechanism on Jupiter. Because the electric 
conductivity increases continuously with depth, it is dif- 
ficult to estimate where in the atmosphere the MHD drag 
acts. For the sake of argument here, we assume that sub- 
stantial MHD drag is confined within ^0.96 Jupiter radii 
(pressures greater than ^10 5 bar). This corresponds to 
the estimated maximum depth up to which zonal flows 
with speeds similar to those in the upper troposphere 
could extend before the energy dissipation by the MHD 
drag would violate the constraint that it cannot exceed 
Jupiter's luminosity (Liu et al. 2008). Precisely at which 
depth the MHD drag acts is not essential for our argu- 
ments, but it is crucial that there is an outer atmospheric 
shell in which the flow is effectively frictionless. 

Because the MHD drag acts at great depth, deep- 
atmosphere dynamics must be considered in linking it 
to dynamics in the upper troposphere. It is well known 
that in a thin atmospheric shell, in the zonal, temporal, 
and vertical mean in a statistically steady state, drag on 
the zonal flow balances any transport of angular momen- 
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FIGURE 2: Schematic of Jupiter's envisaged zonal flow and 
mean meridional circulations. Shown is one quadrant in the 
meridional plane (not to scale). Colors indicate the zonal flow 
(yellow/red for prograde and cyan/blue for retrograde flow). 
Contours indicate the mass flux streamfunction of the mean 
meridional circulation (solid for clockwise rotation and dashed 
for counterclockwise rotation). We have omitted mean merid- 
ional circulations in the equatorial region in which cylinders 
concentric with the planet's spin axis do not intersect the layer 
with MHD drag; they likely have a more complex structure 
than those in higher latitudes. 

turn into or out of a latitude band (e.g., Peixoto and Oort 
1992, chapter 11); that is, the zonally, temporally, and 
vertically averaged angular momentum balance is 

(S e + S m )=r x (D), (13) 

where (■) denotes a mass-weighted vertical mean and T> 
is the drag force per unit mass on the zonal flow, for 
example, T> = —ku with relaxation coefficient k for 
Rayleigh drag. This generalizes to a deep atmosphere 
if r± = r cos <\> is taken to be the actual distance to 
the planet's spin axis, with distance to the planet's cen- 
ter r, and if the vertical mean is understood as a mean 
along surfaces of constant planetary angular momen- 
tum per unit mass Mq = Q,r\. Such angular momen- 
tum surfaces are vertical in the thin-shell approximation 
(r = a = const) but are parallel to the planet's spin 
axis in a deep atmosphere. That is, the zonal and vertical 
mean in a thin atmosphere must be replaced by a mean 
along cylinders concentric with the planet's spin axis in 
a deep atmosphere (see appendix A). 

In the equatorial region in which cylinders concentric 
with the planet's spin axis do not intersect the layer with 
MHD drag, there is effectively (absent other dissipation 



mechanisms except weak viscous dissipation) no cylin- 
drically averaged drag (D) on the zonal flow. If substan- 
tial MHD drag is confined within ^0.96 Jupiter radii, 
the region of vanishing cylindrically averaged MHD 
drag extends in the outer atmosphere from the equator 
to ^16° latitude in each hemisphere and projects from 
there downward along cylinders (see the schematic in 
Fig. 2). In this region, the angular momentum balance 
(13) implies that any convergence of eddy angular mo- 
mentum fluxes in the upper troposphere can only be bal- 
anced by divergence of eddy angular momentum fluxes 
in other layers and/or by divergence of relative angu- 
lar momentum fluxes owing to mean meridional circu- 
lations. The mean meridional circulations likely can- 
not link equatorial convergence of eddy angular momen- 
tum fluxes in the upper troposphere to the MHD drag at 
depth, as considerations of Rossby numbers show. Mean 
flows in the meridional direction can only lead to sub- 
stantial divergence of relative angular momentum fluxes 
where the Rossby number Ro = U/\fL\, with merid- 
ional length scale of zonal-flow variations L, is order 
one or greater. The region where Ro > 1 has meridional 
half-width < L e (footnote 6), even at depths at which 
the zonal flow speed may exceed its upper-tropospheric 
value but remains constrained by the bound (12). It ex- 
tends to < 8° latitude in the outer atmosphere, so the 
latitude band in which Ro > 1 is narrower than the lati- 
tude band of vanishing cylindrically averaged MHD drag 
provided there is no MHD drag in an outer atmospheric 
shell at least ~0.01 Jupiter radii thick. Similarly, mean 
flows in the cylindrically-radial direction (perpendicular 
to the planet's spin axis) can only lead to substantial di- 
vergence of relative angular momentum fluxes where the 
Rossby number Roj^ = U /(2SIL±), with cylindrically- 
radial length scale of zonal-flow variations L±, is or- 
der one or greater (see appendix A). With velocity scale 
U ~ 200 ms -1 and length scale L± ~ 2000 km, corre- 
sponding to the radial length scale in the equatorial plane 
of a jet that extends from the equator to ^13° latitude 
in the outer atmosphere and projects downward along 
cylinders, this Rossby number is Roj^ ~ 0.3; it is even 
smaller in higher latitudes where the zonal-flow velocity 
must be smaller because it cannot substantially exceed 
its upper-tropospheric values of order 10 ms -1 without 
violating the constraint that the energy dissipation by the 
MHD drag cannot exceed Jupiter's luminosity (Liu et al. 
2008). Therefore, mean meridional circulations can re- 
distribute angular momentum only within a region that is 
narrower than that in which cylinders concentric with the 
planet's spin axis do not intersect the layer with MHD 
drag, provided there is no MHD drag in an outer at- 
mospheric shell at least ~0.01 Jupiter radii thick. So 
they likely cannot link equatorial convergence of eddy 
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angular momentum fluxes resulting from meridional ra- 
diation of convectively generated Rossby waves to the 
MHD drag at depth. To the extent that divergence of 
eddy angular momentum fluxes outside the layer of con- 
vective Rossby wave generation and divergence of rela- 
tive angular momentum fluxes owing to mean meridional 
circulations merely compensate the acceleration of the 
cylindrically averaged zonal flow, generation of equato- 
rial superrotation by Rossby wave radiation seems in- 
evitable. 

In higher latitudes, where cylinders concentric with 
the planet's spin axis do intersect the layer with MHD 
drag, mean meridional circulations link dynamics in the 
upper troposphere to the MHD drag at depth. The 
Rossby numbers Ro and Roj^ in higher latitudes are 
small, both in the upper troposphere and in deeper lay- 
ers. Therefore, the relative angular momentum flux di- 
vergence owing to eddies dominates that owing to mean 
meridional circulations, and the local angular momen- 
tum balance (10) in the upper troposphere, where there 
is effectively no drag on the zonal flow, reduces to 

fv x r ± ^S e . (14) 

It follows that in regions of eddy angular momentum flux 
convergence (prograde jets), the mean meridional mass 
flux is equatorward; in regions of divergence (retrograde 
jets), it is poleward. Between the layer with substantial 
eddy angular momentum fluxes and the layer with MHD 
drag, provided the divergence of convective angular mo- 
mentum fluxes (Reynolds stress) is negligible, the mean 
meridional circulations are unaffected by zonal torques 
(S e w r±T> w 0); the local angular momentum balance 
(10) in a form suitable for a deep atmosphere can then 
be expressed as 

u-VM«0, (15) 

where M — Mq,+M u is the angular momentum per unit 
mass. It follows that the mean meridional circulations 
extend downward along surfaces of constant angular mo- 
mentum per unit mass M (e.g., Haynes et al. 1991). 
Approximately, these angular momentum surfaces are 
again cylinders concentric with the spin axis because 
small Rossby numbers mean that the angular momen- 
tum M is dominated by its planetary component Mq. 
Irrespective of the depth at which the MHD drag acts, 
the mean meridional circulations must extend downward 
to and must close where the drag allows mass fluxes to 
cross angular momentum surfaces (Haynes et al. 1991; 
O'Gorman and Schneider 2008). There, the Coriolis 
torque on any mass flux component normal to angular 
momentum surfaces is balanced by the MHD drag on 
the zonal flow (Ekman balance), 

u • VM n « r±T>. (16) 



Figure 2 sketches the resulting mean meridional circula- 
tions. 

Ekman balance (16) at depth requires that the zonal 
flow is prograde {T> < 0) where the mean meridional 
mass flux has a component toward the spin axis and ret- 
rograde (D > 0) where it has a component away from 
the spin axis. Given the correlation between the zonal 
velocity and the convergence of eddy angular momen- 
tum fluxes in the upper branches of the mean meridional 
circulations, this implies that the signs and zeros of the 
zonal flow, like the mean meridional circulations, project 
downward approximately along cylinders, as sketched in 
Fig. 2. Thermal wind balance then constrains the thermal 
structure of the atmosphere below the layer with substan- 
tial eddy angular momentum fluxes [see, e.g., Ingersoll 
and Pollard (1982) and Kaspi (2008) for thermal wind 
equations for deep atmospheres]. The mean meridional 
circulations adjust entropy gradients and the zonal flow 
in lower layers such that they satisfy, in a statistically 
steady state, the constraints that (i) angular momentum 
flux convergence or divergence and the MHD drag on the 
zonal flow balance upon averaging over cylinders, and 
(ii) the zonal flow is in thermal wind balance with the 
entropy gradients (see, e.g., Haynes et al. 1991). These 
dual constraints generally cannot be satisfied, as is often 
assumed in Jupiter models, with entropy gradients that 
vanish throughout the deep atmosphere and with a corre- 
sponding zonal flow without shear in the direction of the 
spin axis (Taylor columns). The zonal flow speed within 
~0.96 Jupiter radii is constrained to be smaller than than 
that in the prograde off-equatorial jets in the upper tro- 
posphere because otherwise the energy dissipation by 
the MHD drag would exceed Jupiter's luminosity (Liu 
et al. 2008). Therefore, the zonal flow generally must 
be sheared and must be associated with nonzero entropy 
gradients. For example, it is conceivable that convection 
deep in Jupiter's atmosphere homogenizes entropy along 
angular momentum surfaces but that there are nonzero 
(albeit possibly weak) entropy gradients normal to them, 
such that the mean thermal structure at depth is neutral 
with respect to slantwise convective, or symmetric, in- 
stability. 

5. General circulation model 

To demonstrate the viability of the mechanisms pro- 
posed, we constructed a GCM of Jupiter's outer atmo- 
sphere. Simulating Jupiter's deep atmosphere in a man- 
ner that is consistent with its measured energy balance is 
computationally prohibitive. Instead, our GCM is based 
on the hydrostatic primitive equations for a dry ideal- 
gas atmosphere in a thin spherical shell with Jupiter's 
radius, rotation rate, gravitational acceleration, and ther- 
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mal properties. The shell extends from the top of the at- 
mosphere to a lower boundary with mean pressure 3 bar. 
Within this pressure range, about 90% of the solar radi- 
ation incident at the top of Jupiter's atmosphere is ab- 
sorbed or scattered back, and latent heat release in phase 
changes of water, ammonia, and hydrogen sulfide is neg- 
ligible: moist adiabatic and dry adiabatic temperature 
lapse rates are nearly indistinguishable (Showman and 
Ingersoll 1998; Ingersoll et al. 2004). Nonetheless, the 
idealization of focusing on a thin shell in Jupiter's outer 
atmosphere means that we are not able to resolve details 
of the coupling between the flow in the outer atmosphere 
and that at depth, beyond the previously discussed con- 
straints on this coupling implied by the angular momen- 
tum balance. 

The GCM uses the spectral-transform method in the 
horizontal (resolution T213) and finite differences in the 
vertical (30 levels). Radiative transfer is represented 
as that in a homogeneous gray atmosphere, with ab- 
sorption and scattering of solar radiation and absorp- 
tion and emission of thermal radiation. Optical prop- 
erties of the atmosphere are so specified that the ideal- 
ized representation of radiative transfer is qualitatively 
consistent with measured radiative fluxes inside Jupiter's 
atmosphere (Sromovsky et al. 1998) and quantitatively 
consistent with the measured energy balance at the top 
of the atmosphere (Hanel et al. 1981). At the GCM's 
lower boundary, a spatially uniform and temporally con- 
stant intrinsic heat flux is imposed. If the intrinsic heat 
flux is sufficiently strong to destabilize the atmospheric 
thermal stratification, the convection that ensues is rep- 
resented by a quasi-equilibrium convection scheme that 
relaxes temperatures in statically unstable parts of atmo- 
spheric columns to a profile with neutral static stability 
(Schneider and Walker 2006). Near the lower boundary, 
Rayleigh drag in the horizontal momentum equations is 
chosen as an idealized linear representation of effects of 
the MHD drag deep in Jupiter's atmosphere on the flow 
in the outer atmosphere. To represent the downward pro- 
jection along cylinders in a deep atmosphere in the thin- 
shell approximation with vertical planetary angular mo- 
mentum surfaces, we use a Rayleigh drag coefficient that 
is constant along the GCM's lower boundary poleward of 
16.3° latitude and rapidly decreases to zero at lower lat- 
itudes, corresponding to the assumption that substantial 
MHD drag is confined within 0.96 Jupiter radii. Above 
the layer with Rayleigh drag, highly scale-selective hori- 
zontal hyperdiffusion at small scales (Smith et al. 2002), 
representing subgrid-scale processes, is the only fric- 
tional process; in particular, there is no vertical viscous 
transfer of momentum or heat. All forcings and bound- 
ary conditions are temporally constant and zonally and 
hemispherically symmetric, so the simulations have sta- 



tionary and zonally and hemispherically symmetric flow 
statistics. See appendix B for details of the GCM and 
simulations. 

6. Jupiter simulation 

a. Upper-tropospheric dynamics 

A simulation with Jupiter's solar constant, with insola- 
tion at the top of the atmosphere varying with the cosine 
of latitude, and with an intrinsic heat flux of 5.7 W m~ 2 
(Gierasch et al. 2000) at the lower boundary reproduces 
the large-scale features of Jupiter's upper-tropospheric 
jets. It produces a prograde equatorial jet and an alternat- 
ing sequence of retrograde and prograde off -equatorial 
jets, with speeds and widths similar to those on Jupiter 
(Fig. la). The jets maintain their speeds and struc- 
ture over time in the statistically steady state of the 
simulation; however, there is some low-frequency vari- 
ability, most pronounced in low latitudes, with decadal 
and possibly longer timescales. The equatorial jet is 
wider and stronger than the off-equatorial jets and re- 
sembles Jupiter's equatorial jet, though it does not ex- 
hibit Jupiter's local velocity minimum at the equator. 
As on Jupiter, retrograde off-equatorial jets are weaker 
than prograde jets, except for the first retrograde jets off 
the equator; in high latitudes, some retrograde jets are 
more manifest as local zonal velocity minima than as 
actual retrograde flow in the upper troposphere. (The 
retrograde jets are stronger at lower levels, but they are 
still weaker than the prograde jets within our GCM do- 
main; see Fig. 4 below.) Prograde off-equatorial jets are 
sharper than retrograde jets, consistent with them being 
barotropically more stable (Rhines 1994). The merid- 
ional gradient of absolute vorticity in the upper tropo- 
sphere is generally positive but, particularly in shorter- 
term averages, reverses in the sharpest retrograde jets, 
where it can approach — f}/2 — similar to, but of smaller 
magnitude than, the reversed absolute vorticity gradi- 
ents in Jupiter's retrograde jets (Ingersoll et al. 1981; 
Read et al. 2006). That Jupiter's retrograde jets vio- 
late sufficient conditions for linear barotropic stability 
in the upper troposphere has been much discussed and is 
sometimes deemed puzzling (e.g., Ingersoll et al. 2004; 
Vasavada and Showman 2005). However, as demon- 
strated here and previously, for example, by Williams 
(2002) and Kaspi and Flierl (2007) for baroclinic flows 
and by Marston et al. (2008) for barotropic flows, it is 
not necessary that statistically steady states of forced- 
dissipative flows satisfy barotropic stability conditions 
for unforced and non-dissipative flows. 

The speeds and widths of the off-equatorial jets in- 
crease with decreasing Rayleigh drag coefficient, sim- 
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ilar to the drag dependence of two-dimensional turbu- 
lent flows described by Smith et al. (2002) and Danilov 
and Gurarie (2002). Because the drag coefficient is 
poorly constrained by data, we empirically chose a value 
(0.05 day^ 1 with lday = 86400 s) that resulted in a 
good fit to Jupiter's off-equatorial jet speeds and widths. 
We also experimented with different widths of the re- 
gion of vanishing Rayleigh drag. Doubling the width 
of this region such that it extends to 33° latitude in 
each hemisphere leads to a prograde equatorial jet that is 
^35 ms -1 stronger at the equator and whose first zero 
off the equator is ^5° farther poleward, giving a closer 
match to Jupiter's equatorial jet; however, it also leads 
to the first retrograde jets off the equator being consider- 
ably stronger and wider than those on Jupiter. Nonethe- 
less, it is clear that the width of the region of vanishing 
drag does not alone control the width of the prograde 
equatorial jet. But having a sufficiently wide region of 
vanishing or very weak drag around the equator is neces- 
sary to obtain equatorial superrotation with intrinsic heat 
fluxes comparable to Jupiter's; if the drag coefficient in 
the equatorial region is similar to that in higher latitudes, 
considerably stronger intrinsic heat fluxes are necessary 
to generate equatorial superrotation (see appendix B for 
details). 

Instantaneous flow fields in the upper troposphere pro- 
vide the first evidence that the equatorial superrotation 
is indeed a consequence of meridional radiation of con- 
vectively generated Rossby waves. Although the im- 
posed intrinsic heat flux at the lower boundary is spa- 
tially uniform, fluctuations in the horizontal divergence 
associated with convective heating fluctuations are pri- 
marily confined to an equatorial latitude band and there 
are typically modulated by large-scale waves; the hor- 
izontal flow away from the equator is geostrophic and 
thus nondivergent to leading order (Fig. 3a). Root-mean- 
square (rms) fluctuations in the horizontal divergence are 
maximal at the equator, where they reach ~1 x 10~ 5 s _1 ; 
they decay to half their maximum value at ^4° latitude, 
consistent with our scaling estimates (footnote 1). In the 
vertical, rms fluctuations in the horizontal divergence are 
maximal in the upper troposphere, below the top of the 
convective outflows (near and slightly above the level at 
which they are shown in Fig. 3). Similarly, the Rossby 
wave source (5) is also primarily confined to an equato- 
rial latitude band (Fig. 3b). In the vertical, rms fluctua- 
tions in the Rossby wave source (5) are maximal in the 
upper troposphere, like the horizontal divergence fluctu- 
ations. In the horizontal, they are maximal ^2° off the 
equator, where they can be generated by vortex stretch- 
ing and reach <~1.5 x 10~ 10 s~ 2 ; they decay to half their 
maximum value at ^7° latitude, consistent with our scal- 
ing estimates (footnote 3). While fluctuations in the hor- 



izontal divergence and in the Rossby wave source (5) 
are primarily confined to the equatorial region, the vor- 
ticity of the horizontal flow shows eddies at all latitudes 
(Fig. 3c). The vorticity fluctuations are of similar mag- 
nitude as divergence fluctuations near the equator but 
are much larger than divergence fluctuations away from 
the equator — consistent with eddy generation by convec- 
tive heating fluctuations near the equator and eddy gen- 
eration by baroclinic instability away from the equator. 
The vorticity field clearly shows the shear zones between 
the zonal jets as well as smaller-scale coherent vortices, 
though no large coherent vortices such as the Great Red 
Spot. It is possible that large vortices such as the Great 
Red Spot form spontaneously but would require longer 
integration times than we can achieve in our simulation, 
or that deep-atmosphere dynamics not captured in our 
simulation are important for their formation and stabil- 
ity. That the zonal jets are present and coherent at every 
instant, not only upon averaging, is most clearly evident 
in the zonal velocity field, which also shows the equato- 
rial waves recognizable in the other flow fields, as well 
as undulations of off-equatorial jets (Fig. 3d). Anima- 
tions (available at www.gps.caltech.edu/~4apio/papers/) 
show that the equatorial waves, organized into large 
wave packets, exhibit retrograde phase velocities, con- 
sistent with them being Rossby waves. The retrograde 
tilt of their phase lines away from the equator (Fig. 3d) 
indicates that they transport angular momentum toward 
the equator (cf. Peixoto and Oort 1992, chapter 11). 

b. Vertical structure and angular momentum fluxes 

The vertical structure of the zonal flow in the simula- 
tion indicates preferential baroclinic eddy generation in 
prograde off-equatorial jets and is consistent with what 
is known about Jupiter's equatorial jet in lower layers 
(Fig. 4a). The speed of the prograde equatorial jet in- 
creases with depth, for example, at the equator, from 
^80 ms -1 between 0.25 and 1 bar to ^160 ms _1 at 
3 bar. This is similar to the increase in zonal flow 
speed seen on Jupiter along the Galileo probe descent 
path at 6.4°N (Atkinson et al. 1998), which lies within 
Jupiter's equatorial jet and so may be comparable to 
latitudes closer to the equator in our simulation with a 
narrower equatorial jet. The speed of the off-equatorial 
zonal flow decreases with depth in prograde jets and in- 
creases, with weaker vertical shear, in retrograde jets. 
For example, the maximum prograde jet speed poleward 
of 20° latitude decreases from 30 ms -1 at 0.65 bar to 
23ms _1 at 3 bar, whereas the maximum retrograde jet 
speed increases from 13 m s _1 at0.65 bartol5ms _1 at 
3 bar. This implies that the prograde off-equatorial jets 
are baroclinically more unstable than the retrograde jets. 
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FIGURE 3: Flow fields at 0.65 bar at one instant in Jupiter simulation, (a) Horizontal divergence, (b) Rossby wave source (5). (c) 
Relative vorticity of horizontal flow, (d) Zonal velocity. The instant shown is within the period for which the mean zonal flow is 
shown in Fig. la. 



The thermal structure of the atmosphere is in thermal 
wind balance with the zonal flow, shows the signature of 
convection penetrating into the upper troposphere, and 
is consistent with what is known about Jupiter's thermal 
structure (Fig. 4b). For example, meridional temperature 
gradients along isobars are poleward (reversed) within 
the equatorial jet, consistent with thermal wind balance 
and zonal flow speeds increasing with depth. Away from 
the equator, they are equatorward in prograde jets but are 
poleward or close to zero in many retrograde jets, con- 
sistent with the opposite signs of the vertical shear in 



prograde and retrograde jets. The atmosphere is stably 
stratified above ~0.5 bar, with a temperature inversion 
above ~0.2 bar and with a level of minimum temper- 
ature lapse rate (maximum Brunt-Vaisala frequency) at 
~0.15 bar that may be identified with the tropopause. 
These features of the thermal structure are qualitatively 
and to a large degree quantitatively consistent with the 
available data for Jupiter's upper troposphere (Simon- 
Miller et al. 2006). A possible quantitative discrepancy 
may be that in the layer between 0.1 and 0.5 bar, which 
contains the thermal emission level on Jupiter as well 
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FIGURE 4: Flow fields in the latitude-pressure plane in Jupiter simulation, (a) Zonal flow (contours) and divergence div(u'v' r±) 
of meridional eddy angular momentum fluxes (colors) . Gray contours for zonal flow speeds between 5 and 20ms -1 , with a contour 
interval of 5ms -1 ; black contours for zonal flow speeds of 40ms" 1 or above, with a contour interval of 20 ms _1 . Solid contours 
for prograde flow and dashed contours for retrograde flow, (b) Temperature (contours, contour interval 10 K) and Brunt- Vaisala 
frequency N (colors). Shown in this and subsequent figures are zonal and temporal means over the same 1500 days for which the 
zonal flow at 0.65 bar is shown in Fig. la. The green part of the latitude axis marks the latitudes with Rayleigh drag. The graphs 
are truncated at 0.1 bar at the top, but the uppermost full level of the GCM has a mean pressure of 0.05 bar. 



as in the simulation, pole-equator temperature contrasts 
along isobars may be a few Kelvin larger in the simu- 
lation than on Jupiter. (If it does not result from data 
and retrieval limitations, the lack of vertical resolution 
and of representation of hazes in our GCM may ac- 
count for this discrepancy.) Beneath ~0.5 bar, the at- 
mosphere in the simulation is nearly neutrally stratified 
by convection. As a consequence of the nearly neutral 
stratification, meridional entropy gradients in the con- 
vective layer are nearly invariant in the vertical. Close 
to the equator, large-scale vertical entropy fluxes lead 
to a weakly stable thermal stratification above ^2 bar, 
consistent with Jupiter's stratification along the Galileo 
probe descent path (Magalhaes et al. 2002). The external 
gravity wave speed implied by the thermal stratification 7 
up to 0.15 bar is ^480 ms -1 at the equator, decreases 

7 We estimated the gravity wave speed as c = N p dp, where 
p s is the pressure at the lower boundary, pt = 0.15 bar is the upper 



to ^-400 ms -1 at 5° latitude, and from there decreases 
gradually but not monotonically to ~300ms _1 at 80° 
latitude — values roughly consistent with the estimated 
gravity wave speed for Jupiter (Ingersoll and Kanamori 
1995). 

Consistent with prograde jets being baroclinically 
more unstable than retrograde jets and with baroclinic 
eddy generation preferentially in prograde jets, eddies in 
the simulation transport angular momentum meridion- 
ally from retrograde into prograde off-equatorial jets in 
layers above ~1 bar (Fig. 4a). Poleward of 20°, the di- 



boundary of the integration, N p = — (p9)~ 1 d p is a static stabil- 
ity measure, and p is density (Schneider and Walker 2006). Because 
the thermal stratification is (nearly) neutral in lower layers, the gravity 
wave speed depends only weakly on the lower boundary of the inte- 
gration. However, it decreases as the upper boundary is lowered (pt is 
increased) within the stably stratified upper troposphere. For example, 
with an upper boundary of pt = 0.2 bar, one obtains a gravity wave 
speed of ~350 m s — 1 at the equator. 
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vergence of meridional eddy angular momentum fluxes 
in retrograde jets reaches up to 170 m 2 s~ 2 and the con- 
vergence in prograde jets up to 470 m 2 s~ 2 , with ex- 
tremal values attained between 0.2 and 0.6 bar. In the 
equatorial region, apparently owing to meridional radi- 
ation of convectively generated Rossby waves, there is 
strong convergence of meridional eddy angular momen- 
tum fluxes (up to 890 m 2 s~ 2 ) between 0.6 and 1 .4 bar — 
below and near the transition from the weakly or neu- 
trally stratified lower layers to the more stably strati- 
fied upper layers. There is weaker divergence (up to 
660 m 2 s~ 2 ) at higher and lower levels, as well as im- 
mediately off the equator between 0.6 and 2 bar. 
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FIGURE 5: Zonally, temporally, and vertically averaged angu- 
lar momentum balance in Jupiter simulation. Torque per unit 
mass owing to divergence of eddy angular momentum fluxes 
(S e ) (blue), to divergence of mean angular momentum fluxes 
(Sm) (green), and to drag on the zonal flow r±(T>) (orange). 

The angular momentum balance is closed in the man- 
ner discussed in section 4. Away from the equator, 
in the region with Rayleigh drag, the Rossby number 
Ro = U/\fL\ is small, and mean meridional circula- 
tions link divergences and convergences of eddy angular 
momentum fluxes in the upper troposphere to the drag 
at the lower boundary. The mean meridional circula- 
tions satisfy the local angular momentum balance (14) in 
the upper troposphere and extend downward along ver- 
tical lines (planetary angular momentum contours in the 
thin-shell approximation); they form thin-shell analogs 
of the circulations sketched in Fig. 2. (In fact, we pro- 
duced Fig. 2 by projecting winds and mean meridional 
circulations from a thin-shell GCM simulation into a 
thick spherical shell.) Although the eddy angular mo- 
mentum fluxes are confined to a relatively thin layer, the 
jets extend to the lower boundary, with the signs and 
zeros of the zonal flow projecting downward approxi- 
mately along vertical lines (Fig. 4a). This allows the 
drag on the zonal flow to balance any convergence or 
divergence of eddy angular momentum fluxes in the ver- 



tical average, such that the vertically averaged angular 
momentum balance in the limit of small Rossby number, 
(S e ) ~ r± (T>), is approximately satisfied (Fig. 5). In 
the region of vanishing Rayleigh drag, the vertically av- 
eraged angular momentum balance (S e +S m ) « is sat- 
isfied by partially compensating divergences and conver- 
gences of meridional eddy angular momentum fluxes in 
different layers (Fig. 4a), as well as vertical eddy fluxes 
and mean meridional circulations redistributing angular 
momentum between layers and between latitude bands 
within the equatorial region (Fig. 5). The small residual 
in the zonally, temporally, and vertically averaged angu- 
lar momentum balance (Fig. 5) is indicative of sampling 
variability and is associated with low-frequency variabil- 
ity in the mean zonal flow. 

The correlation between the convergence of eddy an- 
gular momentum fluxes and the mean flow evident in 
Figs. 4 and 5 implies conversion of eddy kinetic energy 
to mean flow kinetic energy. Above 1 bar, the mean 
conversion rate from eddy to mean flow kinetic energy 
per unit mass 8 is ^2 x 10 _5 Wkg _1 . This conversion 
rate and the eddy angular momentum fluxes themselves 
are of similar magnitude as those inferred for Jupiter's 
upper troposphere by Ingersoll et al. (1981) and Salyk 
et al. (2006). However, unlike what was assumed in 
those studies, the eddy angular momentum fluxes in 
our simulation have a baroclinic structure and are con- 
fined to a relatively thin layer in the upper troposphere. 
As a consequence, the total conversion rate from eddy 
to mean flow kinetic energy in our simulation is only 
0.08 Wm~ 2 , that is, only 0.5% of the energy uptake by 
the atmosphere (14.0 Wm -2 ) — an order of magnitude 
smaller than suggested by Ingersoll et al. (1981) and 
Salyk et al. (2006). It is also only a fraction (19%) of 
the generation rate of eddy kinetic energy by conversion 
from eddy potential energy (0.41 Wm~ 2 ). The remain- 
der of the eddy kinetic energy generation is balanced by 
dissipation (by the Rayleigh drag and by the hyperdiffu- 
sion representing subgrid-scale processes), similarly as 
in Earth's atmosphere (cf. Peixoto and Oort 1992, chap- 
ter 14). 

c. Eddy scales and turbulence characteristics 

Eddies in the simulation exhibit a broad range of length 
scales, with the dominant scales depending on latitude 
and on whether zonal or meridional velocity fluctua- 
tions are considered. The vertically averaged zonal spec- 
trum of the zonal eddy velocity variance u' 2 exhibits 



We calculated energy conversion rates following Lorenz (1955). 
The conversion rate from eddy to mean flow kinetic energy includes 
the (small) contributions owing to vertical eddy fluxes of angular mo- 
mentum and other terms. 
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FIGURE 6: Mass-weighted vertical mean of zonal spectra of 
eddy velocity variances, (a) Zonal velocity variance u' 2 . (b) 
Meridional velocity variance v' 2 . The velocity variance con- 
touring is logarithmic. The magenta line in (b) marks the 
energy-containing zonal wavenumber m e {4>)> defined as the 
first negative moment of the zonal spectrum of the meridional 
eddy velocity variance. 

largest power over a broad range of relatively small 
zonal wavenumbers m < 10 (Fig. 6a). The zonal eddy 
velocity variance is maximal in the equatorial region, 
where it peaks at to — 4 and retains substantial power 
at yet smaller wavenumbers (see the equatorial waves 
and wave packets in the instantaneous zonal velocity 
field in Fig. 3d). In contrast, the analogous spectrum 
of the meridional eddy velocity variance v' 2 has rela- 
tively well defined maxima at larger zonal wavenumbers 
to > 10, except in polar latitudes, where the maxima 
are at smaller wavenumbers (Fig. 6b). At wavenumbers 
to > 10 away from and m > 5 near the poles, zonal and 
meridional eddy velocity variances are of similar magni- 
tude, indicating approximate isotropy of the eddies. At 
smaller wavenumbers, the zonal exceeds the meridional 
eddy velocity variance, indicating that the zonal variance 
is associated with anisotropic, predominantly zonal flow 
structures (Boer and Shepherd 1983; Shepherd 1987a). 
This large-scale zonal variance is likely associated with 
variations in the zonal jets themselves, for example, with 
low-frequency undulations. Therefore, the meridional 



eddy velocity variance is more suitable for defining an 
energy-containing turbulent eddy scale. 

The energy-containing zonal wavenumber m e (<f)), de- 
fined as the first negative moment of the vertically 
averaged zonal spectrum of the meridional eddy ve- 
locity variance, 9 closely follows the maximum of the 
zonal spectrum (Fig. 6b). It increases from ^9 near 
the equator to ~25 near 40° latitude and then de- 
creases toward the poles. The associated zonal length 
scale 2na cos </>/m e decreases from ^50000 km near the 
equator to ~ 18000 km near 40° latitude but then varies 
only weakly, between ~ 10000 km and ^20000 km, in 
higher latitudes. Poleward of ~30° latitude, this zonal 
length scale is similar to the meridional separation scale 
between the off-equatorial jets (Fig. la), which suggests 
that the length scale of approximately isotropic eddies 
controls the jet spacing there. 

The change in the energy-containing zonal wavenum- 
ber and length scale between equatorial and off- 
equatorial regions marks a transition between different 
dynamical regimes: 




Latitude 

FIGURE 7: Zonal wavenumbers in Jupiter simulation. Or- 
ange: energy-containing wavenumber m e . Light blue: Rhines 
wavenumber a cos with L p = 7(3 EKE bt 1 '' 4 //3 1/2 ! 

barotropic eddy kinetic energy per unit mass EKEbt, and em- 
pirical constant 7/3 = 1.6. Green dashed: equatorial Rossby 
wavenumber a cos 4>fL e with L e = (c/(3) 1/2 . Green solid: 
extratropical Rossby wavenumber a cos cj>/L x with L x — 
7ac/|/| and empirical constant 7^ = 1.8. The gravity wave 
speed c entering the Rossby wavenumbers is estimated from 
the thermal stratification between the lower boundary and 
0.15 bar (see footnote 7). 



'That is, given the vertically averaged meridional velocity variance 
spectrum E v (m, </>), the energy-containing wavenumber is denned as 
the integral scale 

,,, (Y, m>1 m- 1 E v (m,4>)\- 1 

This integral scale more closely follows the maximum of the zonal 
spectrum than do other low-order moments of the spectrum. 
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• In the equatorial region, the energy-containing 
wavenumber (~9) is approximately equal to 
the wavenumber corresponding to the equatorial 
Rossby radius L e = (c/(3) 1/2 (Fig. 7). 10 Thus, it 
is approximately equal to the wavenumber of long 
equatorial Rossby waves. Consistently, the zonal 
phase velocity spectrum (Randel and Held 1991) of 
the meridional eddy velocity variance exhibits max- 
imum power at retrograde phase velocities between 
—50 and — 100ms , and the analogous spectrum 
of the meridional eddy flux of angular momentum 
peaks at similar phase velocities. 11 This is in agree- 
ment with theory: long equatorial Rossby waves 
are expected to have retrograde phase velocities of 
about — c/3 ~ — 150ms -1 relative to the mean 
zonal flow (Matsuno 1966; Gill 1982, chapter 11.8). 

• Away from the equator, the energy-containing 
wavenumber is approximately equal to the 
wavenumber corresponding to the extratropical 
Rossby radius L x = J x c/\f\, where we fixed the 
empirical constant at j x = 1.8 (Fig. 7). Thus, it 
scales approximately with the expected wavenum- 
ber of the baroclinically most unstable linear waves. 
In middle and high latitudes, the energy-containing 
wavenumber is also closely approximated by the 
wavenumber corresponding to the Rhines scale 
L/3 = 7 / 3EKEbt 1 ^ 4 //? 1 / 2 , obtained by combining 
the barotropic eddy kinetic energy per unit mass 
EKEbt (the eddy kinetic energy of the vertically 
averaged flow) with j3 (Rhines 1975, 1994); we 
fixed the empirical constant at 7^ = 1.6 as in 
O'Gorman and Schneider (2008). The Rhines 
wavenumber even captures the variations in the 
energy-containing wavenumber between retrograde 
and prograde jets; however, equatorward of ~45° 
latitude, the Rhines wavenumber exceeds the 
energy-containing wavenumber (Fig. 7). 

These results are consistent with convectively gener- 
ated equatorial Rossby waves dominating the meridional 
eddy velocity variance in the equatorial region, and baro- 
clinically generated eddies dominating away from the 
equator. 

10 For an approximately isotropic equatorial Rossby wave with 
meridional decay scale L e , we assume the zonal wavelength is 
2nL e , and hence the zonal wavenumber is 2ira cos </>/(2nL e ) = 
acos<j>/ L e (cf. Matsuno 1966). Other zonal wavenumbers in what 
follows, in which O(l) constants in the relevant length scales L can 
be adjusted in any case, are defined analogously as a cos <p/L, that is, 
without factors of 2n. 

"Long waves with retrograde phase velocities are also 
clearly evident in the zonal velocity animations available at 
www.gps.caltech.edu/~tapio/papers/. 
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FIGURE 8: Eddy kinetic energy spectrum and spectral en- 
ergy fluxes as a function of spherical wavenumber. (a) Spec- 
trum of barotropic eddy kinetic energy (eddy kinetic energy 
of vertically averaged flow), with an n -3 power law in spher- 
ical wavenumber n for comparison, (b) Nonlinear spectral 
flux of total (globally integrated) kinetic energy (solid), de- 
composed into components involving interactions with the 
zonal mean (dash-dotted) and involving only eddy-eddy inter- 
actions (dashed). Positive fluxes indicate downscale and neg- 
ative fluxes upscale energy transfer. The spectrum and spec- 
tral fluxes are calculated following Boer (1983) and Shepherd 
(1987a). 

The coincidence of the energy-containing scale with 
the Rossby radius and the jet separation scale away from 
the equator indicates that the off-equatorial jets form 
without an inverse cascade of barotropic eddy kinetic en- 
ergy, or at least without an inverse energy cascade over 
an extended inertial range. The fact that the energy- 
containing scale is larger than the Rhines scale equator- 
ward of ~45° latitude also points to the absence of an 
inverse energy cascade; the coincidence of the energy- 
containing scale with the Rhines scale at higher latitudes 
is no evidence to the contrary (Schneider and Walker 
2006). Indeed, consistent with the absence of an inverse 
energy cascade, the global barotropic eddy kinetic en- 
ergy spectrum as a function of spherical wavenumber n 
does not show the flattening to an n -5 ' 3 power law at 
large scales that would be expected if an inverse energy 
cascade were occurring (Fig. 8a). The nonlinear spec- 
tral flux of total kinetic energy shows upscale transfer 
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of kinetic energy over a range of wavenumbers between 
n^lO and 100 (Fig. 8b), but it is dominated by interac- 
tions that involve the zonal-mean flow (i.e., the m = 
component, dash-dotted line in Fig. 8b). Interactions 
that involve only eddies (i.e., only zonal wavenumbers 
m > 0, dashed line in Fig. 8b) and that could give rise 
to an inverse energy cascade are much weaker, and they 
do not exhibit an extended wavenumber range in which 
their nonlinear spectral flux is constant, as would be ex- 
pected in an inertial range. It is evident that eddy-mean 
flow interactions, in addition to or in place of nonlinear 
eddy-eddy interactions, are fundamental for the forma- 
tion of the off -equatorial jets (O' Gorman and Schneider 
2007, 2008). This situation is similar to that in Earth's 
atmosphere, in which, likewise, eddy-mean flow inter- 
actions are more important than nonlinear eddy-eddy in- 
teractions for the upscale transfer of kinetic energy, and 
there is no extended inertial range (Shepherd 1987a,b). 

7. Control simulations 

To demonstrate that, in our GCM, differential radiative 
heating is indeed responsible for off-equatorial jets and 
intrinsic convective heat fluxes for equatorial superrota- 
tion, we performed one simulation without differential 
radiative heating and one without intrinsic heat fluxes. In 
the simulation without differential radiative heating, an 
intrinsic heat flux of 5.7 Wm -2 is imposed at the lower 
boundary, as in the Jupiter simulation, but insolation at 
the top of the atmosphere is uniform and equal to the 
global-mean insolation in the Jupiter simulation. In the 
simulation without intrinsic heat fluxes, insolation at the 
top of the atmosphere varies with the cosine of latitude, 
as in the Jupiter simulation, but no heat flux is imposed 
at the lower boundary. Other parameters in these control 
simulations are identical to those in the Jupiter simula- 
tion. 

With the intrinsic heat flux but without differen- 
tial radiative heating, a prograde equatorial jet forms, 
flanked by two retrograde jets; however, there are no 
off-equatorial jets (Fig. lb). With differential radiative 
heating but without intrinsic heat fluxes, off-equatorial 
jets form with similar speeds and widths as in the Jupiter 
simulation; however, the equatorial zonal flow is retro- 
grade (Fig. lb). 

The jets that do form in either case again extend to the 
lower boundary, and the thermal structure is in thermal 
wind balance with the zonal flow. In the simulation with- 
out differential radiative heating (Fig. 9, left column), as 
in the Jupiter simulation, the thermal structure shows the 
signature of convection penetrating into the upper tro- 
posphere, with stable stratification above ^0.5 bar and 
nearly neutral stratification beneath. The speed of the 



prograde equatorial jet again increases with depth, albeit 
more weakly than in the Jupiter simulation (cf. Fig. 4a). 
The speed of the equatorial jet is smaller than in the 
Jupiter simulation, probably because the control simu- 
lation lacks disturbances originating outside the equato- 
rial region, which may trigger convective heating fluc- 
tuations and/or may interact with the equatorial mean 
zonal flow in the Jupiter simulation. The atmosphere 
outside the equatorial region is essentially in radiative- 
convective equilibrium. This illustrates that convection 
does not necessarily generate kinetic energy on large 
scales, although the convective heating in the simulation 
fluctuates at all latitudes (cf. Emanuel et al. 1994). In 
the equatorial region, the weakness of horizontal tem- 
perature gradients implies that convective heating fluc- 
tuations cannot generate substantial temperature fluctua- 
tions but instead, if externally forced, immediately gen- 
erate large-scale horizontal divergence and vorticity fluc- 
tuations and thus large-scale kinetic energy (Sobel et al. 
2001). But outside the equatorial region, convective 
heating fluctuations can generate local temperature fluc- 
tuations, which can decay through radiative processes 
and dispersion by small-scale flows that eventually dis- 
sipate, without generating substantial large-scale kinetic 
energy. The assumption commonly made in shallow- 
water models of Jupiter's upper atmosphere — that con- 
vective heating fluctuations even away from the equator 
directly generate large-scale kinetic energy — is in need 
of justification. 

In the simulation without intrinsic heat fluxes (Fig. 9, 
right column), most of the atmosphere is stably strati- 
fied, except for a nearly neutral layer around ^0.8 bar 
in high latitudes. The flow throughout the equatorial re- 
gion is retrograde. The speed of the off-equatorial zonal 
flow decreases with depth in prograde jets and increases, 
with weaker vertical shear, in retrograde jets, similarly 
as in the Jupiter simulation (cf. Fig. 4a). The prograde 
vertical shear of the zonal flow and, with it, equator- 
ward temperature gradients along isobars are generally 
larger in the simulation without intrinsic heat fluxes, ap- 
parently because the baroclinicity of the stably strati- 
fied atmosphere in the control simulation is smaller than 
that of the atmosphere in the Jupiter simulation with a 
nearly neutrally stratified layer; the smaller baroclinic- 
ity implies weaker meridional eddy entropy fluxes and 
thus larger meridional temperature gradients along iso- 
bars (cf. Schneider and Walker 2008). Consistent with 
preferential eddy generation by baroclinic instability in 
prograde jets, there is again meridional eddy angular 
momentum transport from retrograde into prograde off- 
equatorial jets in layers above ~1 bar. 

The control simulations confirm that, in our GCM, 
both differential radiative heating and intrinsic convec- 
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FIGURE 9: Flow fields in the latitude-pressure plane in control simulations. Top row: Zonal flow (contours) and divergence 
div(u'v' r±) of meridional eddy angular momentum fluxes (colors). Bottom row: Temperature (contours) and Brunt- Vaisala 
frequency N (colors). Contour intervals and plotting conventions as in Fig. 4, except that the contour interval for the angular 
momentum flux divergence is halved. Left column: Simulation with intrinsic heat fluxes but with uniform insolation at the top 
of the atmosphere. Right column: Simulation with differential insolation but without intrinsic heat fluxes. Shown are zonal and 
temporal means over the same 900 simulated days for which the zonal flows at 0.65 bar are shown in Fig. lb. 



tive heat fluxes are necessary to reproduce Jupiter's ob- 
served jets and thermal structure. The imposed intrin- 
sic heat flux needs to exceed a threshold value so that 
convection penetrates into the upper troposphere. Cal- 
culations of radiative equilibrium states show that, with 
the representation of radiative processes in our GCM, the 
intrinsic heat flux needs to exceed ~2Wm -2 for the 
troposphere to become statically unstable. Indeed, the 
zonal flow in a simulation in which an intrinsic heat flux 
of 2 Wm~ 2 is imposed is similar to that in the control 
simulation without intrinsic heat fluxes. 

8. Conclusions and implications 

Based on the theory and simulations and consistent with 
the available observational data, we propose that baro- 
clinic eddies generated by differential radiative heat- 
ing are responsible for Jupiter's off-equatorial jets, and 
that Rossby waves generated by intrinsic convective heat 
fluxes are responsible for the equatorial superrotation. 
Mean meridional circulations adjust entropy gradients 
and the zonal flow in lower layers of Jupiter's atmo- 
sphere, such that the zonal flow is in thermal wind bal- 
ance with the entropy gradients, and convergence or di- 
vergence of angular momentum fluxes in the upper tro- 



posphere and the MHD drag on the zonal flow at depth 
balance upon averaging over cylinders concentric with 
the planet's spin axis. As demonstrated by the Jupiter 
simulation, the resulting view of how the zonal flow and 
general circulation are generated and maintained is con- 
sistent with observed large-scale features of Jupiter's jets 
and thermal structure, such as the zonal flow and merid- 
ional temperature variations in the upper troposphere and 
the thermal stratification of the upper troposphere and 
layers beneath. It is also consistent with the observed 
eddy angular momentum fluxes and with energetic con- 
straints indicating that these fluxes are confined to a rel- 
atively thin atmospheric layer. As demonstrated by the 
control simulations, differential radiative heating alone 
can account for the off -equatorial jets, and intrinsic con- 
vective heat fluxes can account for the prograde equato- 
rial jet. However, intrinsic convective heat fluxes alone 
do not necessarily lead to formation of off-equatorial 
jets. 

The theory and simulations predict aspects of the gen- 
eral circulation that have not been observed but that are 
or soon will be observable. For example, the transition 
in energy-containing eddy scale between the equatorial 
region and regions away from the equator (Fig. 7), point- 
ing to different mechanisms of eddy generation, should 
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be observable by tracking cloud features. And we pre- 
dict that the measurements of NASA's upcoming Juno 
mission to Jupiter will be consistent with zonal jets that 
extend deeply into the atmosphere at all latitudes, away 
from the equator up to depths at which the MHD drag 
acts. As already observed in the upper troposphere, we 
expect that also at lower layers the speeds of the pro- 
grade off-equatorial jets decrease with depth, and that 
there are associated equatorward temperature gradients 
along isobars: qualitatively as in Fig. 4a and b, but likely 
not quantitatively so because the jets are expected to ex- 
tend to much greater depth than in our simulations and 
thus likely have weaker shear. Depending on the strength 
of the MHD drag and on the depth at which it acts, the 
speeds of the retrograde off-equatorial jets may decrease 
or increase with depth, with weaker shear than the pro- 
grade jets, and with associated poleward (reversed) or 
weaker equatorward temperature gradients. If the shear 
of the zonal flow can be inferred from measurements, it 
can be used together with observations of the eddy an- 
gular momentum transport and with the implied transfer 
of kinetic energy from eddies to the mean flow to con- 
strain the strength and depth of the MHD drag and thus 
to elucidate dynamics of the deep atmosphere that are 
not amenable to direct measurement. 

The proposed mechanisms are generic and likely act 
in the atmospheres of all giant planets. They suggest, for 
example, that the reason that Saturn's prograde equato- 
rial jet is wider and stronger than Jupiter's may be that 
Saturn's tropospheric gravity wave speed and equatorial 
Rossby radius are greater. The greater depth at which 
MHD drag on Saturn is estimated to act implies a wider 
region over which there is effectively no drag on the 
zonal flow (Liu et al. 2008), thus making a wider equato- 
rial jet possible. The proposed mechanisms also suggest 
that the reason Uranus and Neptune do not exhibit equa- 
torial superrotation may be that their intrinsic heat fluxes 
are not sufficiently strong to lead to convection penetrat- 
ing into the upper troposphere. 
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APPENDIX A 

Average Angular Momentum Balance and Rossby 
Numbers 

An average angular momentum balance follows from 
the balance equation 

dt(pM) + V • (puM) = -d xP + r ±P V (Al) 

for the angular momentum per unit mass about the 
planet's spin axis, M, with longitude (azimuth) A, three- 
dimensional velocity vector u, distance to the planet's 
spin axis r±, and drag force per unit mass on the zonal 
flow, V. The angular momentum per unit mass M — 
Mq+M u consists of a planetary component Mq — flrj_ 
owing to the planetary rotation and a relative component 
M u = ur± owing to the relative zonal velocity u of the 
atmosphere (e.g., Peixoto and Oort 1992, chapter 1 1). In 
this form, the angular momentum balance (Al) is exact 
(up to neglected viscous stresses) and holds irrespective 
of the constitutive laws of the atmosphere. Averaging 
it temporally and zonally (over azimuth and at constant 
r ±) gives, in a statistically steady state, 

V • ( P uM u + pu Mq )=r ±P ^D, ( A2) 

where the overbar denotes a zonal and temporal mean. 
This is the general, coordinate-independent form of the 
zonally and temporally averaged angular momentum 
balance (10) in pressure coordinates in section 3. If 
one additionally averages (A2) along surfaces of con- 
stant Mq (i.e., over lines of constant r± in the merid- 
ional plane) and uses that, in a statistically steady state, 
there can be no net mass flux across such fixed surfaces, 
it follows that only the zonally and temporally averaged 
relative angular momentum fluxes puM u across an Mq 
surface contribute to any net convergence or divergence 
of angular momentum fluxes and any net acceleration or 
deceleration of the zonal flow averaged over the Mq sur- 
face. This acceleration or deceleration must be balanced 
by drag somewhere on the Mq surface; that is, 

V-{p~^M~}=r ± {pV}, (A3) 

where {•} denotes the mean along lines of constant rj_ 
in the meridional plane. This is the general, coordinate- 
independent form of the angular momentum balance 
(13) in section 4, where the relative angular momentum 
flux puM u is decomposed into mean and eddy compo- 
nents and the zonal and temporal mean is taken along 
isobars. In the thin-shell approximation, rj_ = acoscj> 
depends only on latitude (f>; the average over Mq sur- 
faces is the usual zonal and vertical average (e.g., Vallis 
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2006). In a deep atmosphere, r± = r cos <j) is the actual 
distance to the planet's spin axis; the average over Mq 
surfaces is an average along cylinders concentric with 
the spin axis. 

The relative magnitude of relative and planetary an- 
gular momentum advection by Eulerian mean flows is 
indicated by Rossby numbers, with the relevant Rossby 
number depending on the direction of the flows consid- 
ered. For meridional flows (i.e., flows along meridians 
on spheres concentric with the planet's center), the ra- 
tio of relative to planetary angular momentum advection 
scales as d y M u /d y Mn ~ U/\fL\ because d y Mn — 
—r±f, where y is the meridional coordinate (arc length) 
and L <C r± is the meridional length scale of zonal-flow 
variations; thus, this ratio scales with the usual Rossby 
number Ro = U/\fL\. For cylindrically-radial flows 
(i.e., flows perpendicular to the planet's spin axis), the 
ratio of relative to planetary angular momentum advec- 
tion scales as d r± M u /d r± Mn <~ U/(2ft,L±) because 
d r± Mn — 2£lr±, where L± <C r± is the cylindrically- 
radial length scale of zonal-flow variations; thus, this ra- 
tio scales with the Rossby number Roj^ = U/ (2QL±). 



APPENDIX B 

Jupiter GCM 

The GCM is based on the Flexible Modeling System 
developed at NOAA's Geophysical Fluid Dynamics Lab- 
oratory. It uses standard Jupiter parameters (Table 1), 
except where noted, and integrates the hydrostatic prim- 
itive equations for a dry ideal-gas atmosphere in a thin 
spherical shell, with a stress-free upper boundary at zero 
pressure and a lower boundary with constant geopoten- 
tial and with a mean pressure of 3 bar. 

a. Discretization and resolution 

The primitive equations in Bourke's (1974) vorticity- 
divergence form are discretized with the spectral trans- 
form method in the horizontal, finite differences in the 
vertical, and with semi-implicit time-differencing (e.g., 
Durran 1999, chapter 7.6). The horizontal spectral res- 
olution of the GCM is T213 (triangular truncation of 
the spherical harmonics expansion at wavenumber 213), 
with 640 x 320 (longitude x latitude) points on the Gaus- 
sian transform grid. The vertical coordinate is a = p/p s 
(pressure p normalized by pressure at lower boundary 
p s ) and is discretized with 30 equally spaced levels. With 
this vertical discretization, the density varies by two or- 
ders of magnitude from the top to the bottom of the do- 
main. 



b. Radiative transfer 

Radiative transfer is represented as that in a homoge- 
neous gray atmosphere, using the two-stream approxi- 
mation. The top-of-atmosphere (TOA) insolation is im- 
posed as a perpetual equinox with no diurnal cycle, 



-TTOA = COS ( 

7T 



(Bl) 



where F = 50.7 Wm 2 is the solar constant (Levine 
et al. 1977). 

The solar optical depth t s is linear in pressure to repre- 
sent scattering and absorption by a well-mixed absorber, 



(B2) 



where r s o is the solar optical depth at pressure p . We 
assume diffuse incidence of solar radiation at TOA. The 
solar flux F for a semi-infinite scattering and absorbing 
atmosphere then is (Petty 2006) 



F = F T oA(l-r co )cxp(-rr s ) 



(B3) 



where 



T = 2VT^b yjl-ui-y , (B4) 
and the Bond albedo r m can be represented as 



(B5) 



Here, u> is the single-scattering albedo and 7 the asym- 
metry factor, which we chose to be a; = 0.8 and 7 = 
0.204 to give Jupiter's Bond albedo of = 0.343 (Sro- 
movsky and Fry 2002; Hanel et al. 1981). For the solar 
optical depth, we chose r s0 = 3 at p = 3 bar to give 
a radiative flux F qualitatively consistent with measure- 
ments along the Galileo probe descent path in Jupiter, 
except in the region of hazes in the upper atmosphere 
(Sromovsky et al. 1998). With these parameters, the so- 
lar radiative flux at the lower boundary of the GCM is 
less than 9% of the incident flux at the top of the atmo- 
sphere. 

The thermal optical depth 77 is quadratic in pressure 
to represent collision-induced absorption, 



n = no 



(B6) 



where r;o is the thermal optical depth at pressure p . We 
chose tio — 80 at p = 3 bar, again to give radiative 
fluxes qualitatively consistent with measurements along 
the Galileo probe descent path and to give a thermal 
emission level (where 77 ~ 1) in the vicinity of ~0.4 bar 
(Ingersoll 1990; Sromovsky et al. 1998). 
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TABLE 1: Parameters in Jupiter GCM 



r oiaiiicLcr, syniDoi 


Value 


KclclcIlCc 


Planetary radius (at 3 bar), a 


69.86 x 10 6 m 


Guillot (1999) 


Planetary angular velocity, il 


1.7587 x 10" 4 s 


Donivan and Carr (1969) 


Gravitational acceleration, g 


26.0ms" 2 


Lodders and Fegley (1998) 


Specific gas constant, R 


3605.38 J kg" 1 K- 1 


Lodders and Fegley (1998) 


Adiabatic exponent, n 


2/7 




Specific heat capacity, c p = R/k 


12619.0 J kg" 1 K- 1 




Solar constant, Fo 


50.7Wm- 2 


Levine et al. (1977) 


Intrinsic heat flux 


5.7Wm- 2 


Gierasch et al. (2000) 


Bond albedo, r x 


0.343 


Hanel et al. (1981) 


Single-scattering albedo, Cu 


0.8 


Sromovsky and Fry (2002) 



At the lower boundary, energy conservation is im- 
posed: at each grid point, the upward thermal radiative 
flux is set equal to the sum of the downward solar and 
thermal radiative fluxes. 

c. Intrinsic heat flux 

A spatially uniform and temporally constant heat flux 
is deposited in the lowest layer of the GCM to mimic 
Jupiter's intrinsic heat flux. We use a heat flux of 
5.7 W m- 2 (Gierasch et al. 2000), except in the control 
simulation without intrinsic heat fluxes. 

d. Convection scheme 

A quasi-equilibrium convection scheme represents (dry) 
convection. It relaxes temperature profiles toward a 
convective profile with adiabatic lapse rate g/c p « 
2.1 K km- 1 whenever an air parcel lifted adiabati- 
cally from the lowest model level has positive convec- 
tive available potential energy (Schneider and Walker 
2006). The convective relaxation time is 6 h, chosen 
as roughly the time it takes a gravity wave with speed 
c 450 ms- 1 to traverse the equatorial Rossby radius 
(c/ (3) 1 / 2 w 9500 km. The convection scheme conserves 
enthalpy integrated over atmospheric columns; it does 
not transport momentum. It can be viewed as a dry limit 
of the Betts-Miller convection scheme (Betts 1986; Betts 
and Miller 1986). 

e. Drag at lower boundary 

In Jupiter's atmosphere, the conductivity of hydrogen in- 
creases with depth, making a continuous transition from 
very low values in the outer atmosphere to a constant 
value reached at the depth at which hydrogen becomes 
metallic, which occurs at ~1.4Mbar or ~0.84 Jupiter 
radii (Nellis et al. 1996; Liu et al. 2008). Where the 
atmosphere is electrically conducting, the interaction of 



the magnetic field with the flow leads to Ohmic dissipa- 
tion and retards the flow (Liu et al. 2008). 

As an idealized representation of effects of this MHD 
drag on the flow in the outer atmosphere, we use 
Rayleigh drag near the GCM's lower boundary in the 
horizontal momentum equations, 

dtV + ~- = -k{<j),a)v. (B7) 

As in Held and Suarez (1994), the Rayleigh drag coeffi- 
cient k(<fi, a) decreases linearly in a from its value fco(</>) 
at the lower boundary at cr = 1 to zero at <7f, = 0.8, 

k((j), cr) = k (<t>) max ^0, ^ ~ ^ . (B8) 

To represent in the thin-shell approximation the down- 
ward projection along cylinders in a deep atmosphere, 
we use a drag coefficient ko(<f)) that is constant (ko = 
0.05 day -1 ) poleward of 16.3° latitude, corresponding 
to 0.96 planetary radii in a projection onto the equatorial 
plane, and that exponentially decreases to zero at lower 
latitudes (Fig. 10). 

0.05 H 1 I 

> 
a 

2. 

c 

0) 

o 
a) 



o H 1 1 — 1 — i 1 1 

-60° 0° 60° 

Latitude 

FIGURE 10: Rayleigh drag coefficient ko((j>) at the GCM's 
lower boundary. 

The kinetic energy dissipated by the Rayleigh drag is 
returned to the flow as heat to conserve energy. 
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/ Subgrid-scale dissipation 

Horizontal hyperdiffusion in the vorticity, divergence, 
and temperature equations acts at all levels and is the 
only frictional process above the layer with Rayleigh 
drag (c < 0.8). The hyperdiffusion is represented by 
an exponential cutoff filter (Smith et al. 2002), with a 
damping timescale of 2 h at the smallest resolved scale 
and with no damping for spherical wavenumbers smaller 
than 100. 



g. Simulations 

The simulations were spun-up from radiative-convective 
equilibrium temperature profiles with no flow, with small 
random perturbations in temperature and vorticity to 
break the axisymmetry of the initial state. All simula- 
tions were first spun-up at T85 horizontal resolution for 
10000 simulated (Earth) days. The end states of the T85 
simulations were used as initial states of T213 simula- 
tions, which were spun up for at least 18000 additional 
days. Over the first 10000 days of the T213 simulations, 
we experimented with Rayleigh drag parameters near the 
lower boundary to obtain a good fit to Jupiter's observed 
zonal flow. The Rayleigh drag parameters were held 
fixed at the values stated above for at least 8000 days of 
spin-up of the T213 simulations, until statistically steady 
states were reached. 

In the statistically steady states, the global-mean out- 
going thermal radiative flux is within < 0.02 Wm -2 of 
the sum of the global-mean absorbed solar radiative flux 
and the imposed intrinsic heat flux. The vertically in- 
tegrated Rayleigh drag on the zonal flow approximately 
matches the vertically integrated total (mean plus eddy) 
angular momentum flux convergence at each latitude 
(Fig. 5). 

The circulation statistics shown are computed from 
states sampled 4 times daily in the statistically steady 
states of the simulations. Statistics for the Jupiter sim- 
ulation are computed from 1500 simulated days; statis- 
tics for the control simulations are computed from 900 
simulated days. The statistics are first computed on 
the GCM's a surfaces, with the appropriate surface 
pressure-weighting of the averages (e.g., Walker and 
Schneider 2006), and are then interpolated to pressure 
surfaces for display purposes. The divergence of eddy 
angular momentum fluxes is computed as the isobaric di- 
vergence of the eddy fluxes interpolated to pressure sur- 
faces; this divergence does not differ significantly from 
the divergence of eddy angular momentum fluxes on a 
surfaces. 



h. Sensitivity to Rayleigh drag coefficient 

If the Rayleigh drag coefficient ko(4>) is taken to be con- 
stant in latitude and equal to the value we use outside 
the equatorial region (0.05 day -1 ), much higher intrin- 
sic heat fluxes are needed to generate equatorial super- 
rotation. In computationally less demanding simulations 
of a planet with four times Earth's radius but Jovian pa- 
rameters otherwise (as in Table 1), we found that the 
intrinsic heat fluxes had to be of order 100 Wm -2 to 
produce equatorial superrotation for a constant drag co- 
efficient of 0.05 day -1 . If a smaller constant drag coef- 
ficient is chosen (< 0.001 day -1 ), equatorial superrota- 
tion can be generated with intrinsic heat fluxes compa- 
rable to 5.7 Wm -2 , but then the off-equatorial jets are 
much wider and stronger than Jupiter's. The drag in the 
equatorial region has to be sufficiently weak to produce 
equatorial superrotation with the observed intrinsic heat 
flux, and the drag outside the equatorial region has to be 
sufficiently strong to reproduce the speeds and widths of 
Jupiter's off -equatorial jets. 
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